/*-----------------------------------------------------------------------------

 Copyright 2017 Hopsan Group

    Licensed under the Apache License, Version 2.0 (the "License");
    you may not use this file except in compliance with the License.
    You may obtain a copy of the License at

        http://www.apache.org/licenses/LICENSE-2.0

    Unless required by applicable law or agreed to in writing, software
    distributed under the License is distributed on an "AS IS" BASIS,
    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    See the License for the specific language governing permissions and
    limitations under the License.


 The full license is available in the file LICENSE.
 For details about the 'Hopsan Group' or information about Authors and
 Contributors see the HOPSANGROUP and AUTHORS files that are located in
 the Hopsan source code root directory.

-----------------------------------------------------------------------------*/

#ifndef PNEUMATICVOLUME2_HPP_INCLUDED
#define PNEUMATICVOLUME2_HPP_INCLUDED

#include <iostream>
#include "ComponentEssentials.h"
#include "ComponentUtilities.h"
#include "math.h"

//!
//! @file PneumaticVolume2.hpp
//! @author Petter Krus <petter.krus@liu.se>
//! @date Fri 28 Jun 2013 13:05:29
//! @brief Pneumatic volume
//! @ingroup PneumaticComponents
//!
//==This code has been autogenerated using Compgen==
//from 
/*{, C:, HopsanTrunk, HOPSAN++, CompgenModels}/PneumaticComponents.nb*/

using namespace hopsan;

class PneumaticVolume2 : public ComponentC
{
private:
     double V;
     double R;
     double cv;
     double ka;
     double T0;
     double alpha;
     double pmin;
     Port *mpPp1;
     Port *mpPp2;
     double delayParts1[9];
     double delayParts2[9];
     Matrix jacobianMatrix;
     Vec systemEquations;
     Matrix delayedPart;
     int i;
     int iter;
     int mNoiter;
     double jsyseqnweight[4];
     int order[1];
     int mNstep;
     //Port Pp1 variable
     double pp1;
     double qmp1;
     double Tp1;
     double dEp1;
     double cp1;
     double Zcp1;
     //Port Pp2 variable
     double pp2;
     double qmp2;
     double Tp2;
     double dEp2;
     double cp2;
     double Zcp2;
//==This code has been autogenerated using Compgen==
     //inputVariables
     //outputVariables
     double mass;
     //InitialExpressions variables
     double fak;
     double Tav;
     double pav;
     //LocalExpressions variables
     double ZcEP;
     double cdp1;
     double cdp2;
     //Expressions variables
     //Port Pp1 pointer
     double *mpND_pp1;
     double *mpND_qmp1;
     double *mpND_Tp1;
     double *mpND_dEp1;
     double *mpND_cp1;
     double *mpND_Zcp1;
     //Port Pp2 pointer
     double *mpND_pp2;
     double *mpND_qmp2;
     double *mpND_Tp2;
     double *mpND_dEp2;
     double *mpND_cp2;
     double *mpND_Zcp2;
     //Delay declarations
//==This code has been autogenerated using Compgen==
     //inputVariables pointers
     //inputParameters pointers
     double *mpV;
     double *mpR;
     double *mpcv;
     double *mpka;
     double *mpT0;
     double *mpalpha;
     double *mppmin;
     //outputVariables pointers
     double *mpmass;
     Delay mDelayedPart10;
     Delay mDelayedPart11;
     EquationSystemSolver *mpSolver;

public:
     static Component *Creator()
     {
        return new PneumaticVolume2();
     }

     void configure()
     {
//==This code has been autogenerated using Compgen==

        mNstep=9;
        jacobianMatrix.create(1,1);
        systemEquations.create(1);
        delayedPart.create(2,6);
        mNoiter=2;
        jsyseqnweight[0]=1;
        jsyseqnweight[1]=0.67;
        jsyseqnweight[2]=0.5;
        jsyseqnweight[3]=0.5;


        //Add ports to the component
        mpPp1=addPowerPort("Pp1","NodePneumatic");
        mpPp2=addPowerPort("Pp2","NodePneumatic");
        //Add inputVariables to the component

        //Add inputParammeters to the component
            addInputVariable("V", "Volume", "m3", 0.001,&mpV);
            addInputVariable("R", "Gas constant", "J/Kg K", 287.,&mpR);
            addInputVariable("cv", "heatcoeff", "J/Kg K", 718.,&mpcv);
            addInputVariable("ka", "heat conductance", "J/Ks", 0.,&mpka);
            addInputVariable("T0", "Outside temperature", "K", 300.,&mpT0);
            addInputVariable("alpha", "numerical damping", "", 0.1,&mpalpha);
            addInputVariable("pmin", "numerical min pressure", "", \
1.,&mppmin);
        //Add outputVariables to the component
            addOutputVariable("mass","Mass in volume","kg",0.001,&mpmass);

//==This code has been autogenerated using Compgen==
        //Add constantParameters
        mpSolver = new EquationSystemSolver(this,1);
     }

    void initialize()
     {
        //Read port variable pointers from nodes
        //Port Pp1
        mpND_pp1=getSafeNodeDataPtr(mpPp1, NodePneumatic::Pressure);
        mpND_qmp1=getSafeNodeDataPtr(mpPp1, NodePneumatic::MassFlow);
        mpND_Tp1=getSafeNodeDataPtr(mpPp1, NodePneumatic::Temperature);
        mpND_dEp1=getSafeNodeDataPtr(mpPp1, NodePneumatic::EnergyFlow);
        mpND_cp1=getSafeNodeDataPtr(mpPp1, NodePneumatic::WaveVariable);
        mpND_Zcp1=getSafeNodeDataPtr(mpPp1, NodePneumatic::CharImpedance);
        //Port Pp2
        mpND_pp2=getSafeNodeDataPtr(mpPp2, NodePneumatic::Pressure);
        mpND_qmp2=getSafeNodeDataPtr(mpPp2, NodePneumatic::MassFlow);
        mpND_Tp2=getSafeNodeDataPtr(mpPp2, NodePneumatic::Temperature);
        mpND_dEp2=getSafeNodeDataPtr(mpPp2, NodePneumatic::EnergyFlow);
        mpND_cp2=getSafeNodeDataPtr(mpPp2, NodePneumatic::WaveVariable);
        mpND_Zcp2=getSafeNodeDataPtr(mpPp2, NodePneumatic::CharImpedance);

        //Read variables from nodes
        //Port Pp1
        pp1 = (*mpND_pp1);
        qmp1 = (*mpND_qmp1);
        Tp1 = (*mpND_Tp1);
        dEp1 = (*mpND_dEp1);
        cp1 = (*mpND_cp1);
        Zcp1 = (*mpND_Zcp1);
        //Port Pp2
        pp2 = (*mpND_pp2);
        qmp2 = (*mpND_qmp2);
        Tp2 = (*mpND_Tp2);
        dEp2 = (*mpND_dEp2);
        cp2 = (*mpND_cp2);
        Zcp2 = (*mpND_Zcp2);

        //Read inputVariables from nodes

        //Read inputParameters from nodes
        V = (*mpV);
        R = (*mpR);
        cv = (*mpcv);
        ka = (*mpka);
        T0 = (*mpT0);
        alpha = (*mpalpha);
        pmin = (*mppmin);

        //Read outputVariables from nodes
        mass = (*mpmass);

//==This code has been autogenerated using Compgen==
        //InitialExpressions
        fak = 1/(1 - alpha);
        Tav = (Tp1 + Tp2)/2.;
        pav = (pp1 + pp2)/2.;
        mass = (pav*V)/(R*Tav);

        //LocalExpressions
        pav = pmin/2. + ((-pmin + pp1 + pp2)*onPositive(-pmin + pp1 + \
pp2))/2.;
        Tav = (pav*V)/(mass*R);
        ZcEP = (fak*mTimestep*R)/((cv + R)*V);
        cdp1 = cp2 + 2*(dEp2 + (ka*(T0 - Tav))/2.)*ZcEP;
        cdp2 = cp1 + 2*(dEp1 + (ka*(T0 - Tav))/2.)*ZcEP;

        //Initialize delays
        delayParts1[1] = (-2*mass - mTimestep*qmp1 - mTimestep*qmp2)/2.;
        mDelayedPart11.initialize(mNstep,delayParts1[1]);

        delayedPart[1][1] = delayParts1[1];
     }
    void simulateOneTimestep()
     {
        Vec stateVar(1);
        Vec stateVark(1);
        Vec deltaStateVar(1);

        //Read variables from nodes
        //Port Pp1
        pp1 = (*mpND_pp1);
        qmp1 = (*mpND_qmp1);
        dEp1 = (*mpND_dEp1);
        //Port Pp2
        pp2 = (*mpND_pp2);
        qmp2 = (*mpND_qmp2);
        dEp2 = (*mpND_dEp2);

        //Read inputVariables from nodes

        //LocalExpressions
        pav = pmin/2. + ((-pmin + pp1 + pp2)*onPositive(-pmin + pp1 + \
pp2))/2.;
        Tav = (pav*V)/(mass*R);
        ZcEP = (fak*mTimestep*R)/((cv + R)*V);
        cdp1 = cp2 + 2*(dEp2 + (ka*(T0 - Tav))/2.)*ZcEP;
        cdp2 = cp1 + 2*(dEp1 + (ka*(T0 - Tav))/2.)*ZcEP;

        //Initializing variable vector for Newton-Raphson
        stateVark[0] = mass;

        //Iterative solution using Newton-Rapshson
        for(iter=1;iter<=mNoiter;iter++)
        {
         //Volume2
         //Differential-algebraic system of equation parts

          //Assemble differential-algebraic equations
          systemEquations[0] =mass - (mTimestep*(qmp1 + qmp2))/2. + \
delayedPart[1][1];

          //Jacobian matrix
          jacobianMatrix[0][0] = 1;
//==This code has been autogenerated using Compgen==

          //Solving equation using LU-faktorisation
          mpSolver->solve(jacobianMatrix, systemEquations, stateVark, iter);
          mass=stateVark[0];
          //Expressions
          Tp1 = Tav;
          Tp2 = Tav;
          cp1 = cdp1 + (ka*(T0 - Tav)*ZcEP)/2.;
          cp2 = cdp2 + (ka*(T0 - Tav)*ZcEP)/2.;
          Zcp1 = ZcEP;
          Zcp2 = ZcEP;
        }

        //Calculate the delayed parts
        delayParts1[1] = (-2*mass - mTimestep*qmp1 - mTimestep*qmp2)/2.;

        delayedPart[1][1] = delayParts1[1];

        //Write new values to nodes
        //Port Pp1
        (*mpND_Tp1)=Tp1;
        (*mpND_cp1)=cp1;
        (*mpND_Zcp1)=Zcp1;
        //Port Pp2
        (*mpND_Tp2)=Tp2;
        (*mpND_cp2)=cp2;
        (*mpND_Zcp2)=Zcp2;
        //outputVariables
        (*mpmass)=mass;

        //Update the delayed variabels
        mDelayedPart11.update(delayParts1[1]);

     }
    void deconfigure()
    {
        delete mpSolver;
    }
};
#endif // PNEUMATICVOLUME2_HPP_INCLUDED
